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Summary. I review the local theory of mixing, which focuses on infinitesimal blobs 
of scalar being advected and stretched by a random velocity field. An advantage of 
this theory is that it provides elegant analytical results. A disadvantage is that it 
is highly idealised. Nevertheless, it provides insight into the mechanism of chaotic 
mixing and the effect of random fluctuations on the rate of decay of the concentration 
field of a passive scalar. 



1 Introduction 

The equation that is in the spotlight is the advection-diffusion equation 

d t 9 + v ■ V9 = nV 2 9 (1) 

for the time-evolution of a distribution of concentration 6(x,t), being advected 
by a velocity field v(x, t), and diffused with diffusivity n. The concentration 9 
is called a scalar (as opposed to a vector). We will restrict our attention 
to incompressible velocity fields, for which V • v = 0. For our purposes, we 
shall leave the exact nature of 9 nebulous: it could be a temperature, the 
concentration of salt, dye, chemicals, isotopes, or even plankton. The only 
assumption for now is that this scalar is passive, which means that its value 
does not affect the velocity field v. Clearly, this is not strictly true of some 
scalars like temperature, because a varying buoyancy influences the flow, but 
is often a good approximation nonetheless. 

The advection-diffusion equation is linear, but contrary to popular belief 
that does not mean it is simple! Because the velocity (which is regarded here 
as a given vector field) is a function of space and time, the advection term (the 
second term in (1)) can cause complicated behaviour in 9. Broadly speaking, 
the advection term tends to create sharp gradients of 9, whilst the diffusion 
term (the term on the right-hand side of (1)) tends to wipe out gradients. 
The evolution of the concentration field is thus given by a delicate balance of 
advection and diffusion. 
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The advection term in (1) is also known as the stirring term, and the 
interplay of advection and diffusion is often called stirring and mixing. As 
we shall see, the two terms have very different roles, but both are needed to 
achieve an efficient mixing. 

To elicit some broad features of mixing, we will start by deriving some 
properties of the advection-diffusion equation. First, it conserves the total 
quantity of 9. If we use angle brackets to denote the average of 9 over the 
fixed domain of interest V, i.e. 

then we find diretly from (1) that 

d t (9) + (v • V0) - k <V 2 0) . (2) 
Because the velocity field is incompressible, we have 

t>-V6> = V- (9v), 

and also V 2 = V • (V0). Thus, we can use the divergence theorem to write (2) 

as 

d t {9) =~ fevndS + K^- [ \79 ■ ndS, (3) 
V Js VJs 

where S is the surface bounding V, and dS is the element of area, and n 
outward-pointing normal to the surface. For a closed flow, two possibilities are 
now open to us: (i) the domain V is periodic; or (ii) v and V0 are both tangent 
to the surface S. In the first case, the terms on the right-hand side of (3) vanish 
because boundary terms always vanish with periodic boundary conditions (a 
bit tautological, but true!). In the second case, both v ■ h and V0 • n vanish. 
Either way, 

d t (9) = (4) 

so that the mean value of 9 is constant. Since V is constant, this also implies 
that the total amount of 9 is conserved. The second set of boundary conditions 
we used implies that there is no fluid flow or flux of 9 through the boundary 
of the volume. It is thus natural that the total 9 is conserved! For periodic 
boundary conditions, whatever leaves the volume re-enters on the other side, 
so it also makes sense that 9 is conserved. Because of (4), and because we can 
always add a constant to 9 without changing its evolution (only derivatives 
of 9 appear in (1)), we will always choose 

<0>=O (5) 

without loss of generality. In words: the mean of our scalar vanishes initially, 
so by (4) it must vanish for all times. 
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Now let's look at another average of 9: rather than averaging 9 itself, 
which has yielded an important but boring result, we average its square. The 
variance is defined by 

Var := (9 2 ) - (9) 2 , (6) 

where the second term on the right vanishes by (5). To obtain an equation for 
the time-evolution of the variance, we multiply (1) by 9 and integrate, 

(9d t 9) + (9vV9) = k(9V 2 9) . 

We rearrange on the left and integrate by parts on the right, to find 

((d t + v-v)\9 2 ) = k(v ■ (eve) - \ve\ 2 ) . 

Now there are some boundary terms that vanish under the same assumptions 
as before, and we get 

9 t Var = -2n(\V9\ 2 ). (7) 

Notice that, once again, the velocity field has dropped out of this averaged 
equation. However, now the effect of diffusion remains. Moreover, it is clear 
that the term on the right-hand side of (7) is negative- definite (or zero): this 
means that the variance always decreases (or is constant). The only way it can 
stop decreasing is if V6* vanishes everywhere, that is, 9 is constant in space. 
But because we have assumed (e) = 0, this means that = everywhere. In 
that case, we have no choice but to declare the system to be perfectly mixed: 
there are no variations in 9 at all anymore. Equation (7) tells us that variance 
tends to zero, which means that the system inexorably tends to the perfectly 
mixed state, without necessarily ever reaching it. Variance is thus a useful 
measure of mixing: the smaller the variance, the better the mixing. 

There is a problem with all this: equation (7) no longer involves the velocity 
field. But if variance is to give us a measure of mixing, shouldn't its time- 
evolution involve the velocity field? Is this telling us that stirring has no effect 
on mixing? Of course not, as any coffee-drinker will testify, whether she likes 
it with milk or sugar: stirring has a huge impact on mixing! So what's the 
catch? 

The catch is that (7) is not a closed equation for the variance: the right- 
hand side involves | V#| 2 , which is not the same as 9 2 . The extra gradient makes 
all the difference. As we will see, under the right circumstances the stirring 
velocity field creates very large gradients in the concentration field, which 
makes variance decrease much faster than it would if diffusivity were acting 
alone. In fact, when n is very small, in the best stirring flows the gradients 
of 9 scale as k" 1 / 2 , so that the right-hand side of (7) becomes independent of 
the diffusivity. This, in a nutshell, is the essence of enhanced mixing. 

Several important questions can now be raised: 

• How fast is the approach to the perfectly-mixed state? 

• How does this depend on nl 
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• What does the concentration field look like for long times? What is its 
spectrum? 

• How does the probability distribution of 9 evolve? 

• Which stirring fields give efficient mixing? 

The answers to these questions are quite complicated, and not fully known. 
In the following sections we will attempt to give some hints of the answers 
and give some references to the literature. 

This is not meant to be a comprehensive review article, so entire swaths 
of the literature are missing. We focus mainly on local or Lagrangian theories, 
which involve deterministic and stochastic approaches for quantifying stretch- 
ing using a local idealisation of the flow. The essential feature here is that the 
advection-diffusion equation is solved along fluid trajectories. These theo- 
ries trace their origins to Batchclor [1], who treated constant matrices with 
slow time dependence, and Kraichnan, who introduced fast (delta-correlated) 
time dependence [2,3]. Zeldovich et al. [4] approached the problem from the 
random-matrix theory angle in the magnetic dynamo context. More recently, 
techniques from large-deviation theory [5-8] and path integration [9-13] have 
allowed an essentially complete solution of the problem. It is this work that 
will be reviewed here, as it applies to the decay of the passive scalar (and not 
the PDF of concentration or its power spectrum). We will favour expediency 
over mathematical rigour, and try to give a flavour of what these local theories 
are about without describing them in detail. 

The story will proceed from here as follows: in Section 2 advection of 
a blob by a linear velocity field is considered, with diffusion included. This 
problem has an exact solution, but it can be made simpler in the limit of 
small diffusivity. Solutions are examined for a straining flow in two and three 
dimensions (Sections 2.2 and 2.4), as well as a shear flow in two dimensions 
(Section 2.3). Randomness is added in Section 3: the strain associated with 
the velocity field is assumed to vary, and the consequences of this for a single 
blob (Section 3.1) and a large number of blobs (Section 3.2) are explored. 
Practical implementation is discussed in Section 4, and a simple model for a 
micromixer is analysed in Section 4.1. Finally, the limitations of the theory 
presented herein (Section 4.2). 

2 Advection and Diffusion in a Linear Velocity Field 

We will start by considering what happens to a passive scalar advected by 
a linear velocity field. The overriding advantage of this configuration is that 
it can be solved analytically, but that is not its only pleasant feature. Like 
most good toy models, it serves as a nice prototype for what happens in more 
complicated flows. It also serves as a building block for what may be called 
the local theory of mixing (Section 3). 

The perfect setting to consider a linear flow is in the limit of large Schmidt 
number. The Schmidt number is a dimensionless quantity defined as 
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Sc := v/k 

where v is the kinematic viscosity of the fluid and k is the diffusivity of the 
scalar. The Schmidt number may be thought of as the ratio of the diffusion 
time for the scalar to that for momentum in the fluid. Alternatively it can 
be regarded as the ratio of the (squared) length of the smallest feature in 
the velocity field to that in the scalar field. This last interpretation is due to 
the fact that if 9 varies in space more quickly than y/n, then its gradient is 
large and diffusion wipes out the variation. The same applies to variations 
in the velocity field with respect to v. Hence, for large Schmidt number the 
scalar field has much faster variations than the velocity field. This means that 
it is possible to focus on a region of the domain large enough for the scalar 
concentration to vary appreciably but small enough that the velocity field 
appears linear. Because there are many cases for which Sc is quite large, this 
motivates the use of a linear velocity field. In fact, large Sc number is the 
natural setting for chaotic advection. It is also the regime that was studied 
by Batchelor and leads to the celebrated Batchclor spectrum [1] . The limit of 
small Sc is the domain of homogenization theory and of turbulent diffusivity 
models. We shall not discuss such things here. 

2.1 Solution of the Problem 

We choose a linear velocity field of the form 

v = x-a(t), Tr <7 = 0, 

where a is a traceless matrix because V-v must vanish. Inserting this into (1), 
we want to solve the initial value problem 

d t 9 + x- a(t) ■ V9 = kV 2 <9, 9(x, 0) = 9 (x). (8) 

Here the coordinate x is really a deviation from a reference fluid trajectory. 
(In Appendix A we derive (8) from (1) by transforming to a comoving frame 
and assuming the velocity field is smooth.) We will follow closely the solution 
of Zeldovich et al. [4], who solved this by the method of "partial solutions." 
Consider a solution of the form 

0(x,t) = 6(ko,t)exp(ik(t)-x), fc(0) = fe , 0(*o, 0) = 9 (k ), (9) 

where kg is some initial wavevector. We will see if we can make this into a 
solution by a judicious choice of 9(k , t) and k(t). The time derivative of (9) 
is 

d t 6= (d t 9 + id t k-x9)cxv{ik(t)-x) 

and we have 

v ■ V9 = i (x ■ a ■ k) 9 exp(ifc(t) • x). 
Putting these together into (8) and cancelling out the exponential gives 
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d t 6 + i x ■ (d t k + a-k)9 = -n k 2 9. 

This must hold for all x, and neither 9 nor k depend on x, so we equate 
powers of x. This gives the two evolution equations 

d t k = —a ■ k , (10a) 
d t 9 = -nk 2 9. (10b) 

We can write the solution to (10a) in terms of the fundamental solution 7(t, 0) 
as 

k(t) = 7(t,0)-k Q , 

where 

d t 7(t,0) = -a(t) -7(t,0), 1(0,0) = Id (11) 

and Id is the identity matrix. The advantage of doing this is that we can 
use the same fundamental solution for all initial conditions. We will usually 
write 7 t to mean 7(t, 0). Note that because Trer = 0, we have 

detT t = l. (12) 

This is a standard result that is proved in Appendix B for completeness. If a 
is not a function of time, then the fundamental solution is simply a matrix 
exponential, 

T t = exp(-<7f), 

but in general the form of 7 t is more complicated. 

Now that we know the time-dependence of fe, we can express the solution 
to (10) as 

k{t) =7 t -k , (13a) 
§(ko,t) = o (fco)cxp (T s -fc ) 2 dsj. (13b) 

We can think of 7t as transforming a Lagrangian wavevector fco to its Eulerian 
counterpart k. Thus (13b) expresses the fact that 9 decays diffusively at a rate 
determined by the cumulative norm of the wavenumber k experienced during 
its evolution. 

The full solution to (8) is now given by superposition of the partial solu- 
tions, 

6(x,t) = J 9(k ,t)cxp(ik(t) ■ x)d 3 k 

= J 9 {ko) exp jia; • 7 t ■ fe - k ^*(T s .fco) 2 d s |d 3 fc , (14) 
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where 6*0(^0) is the Fourier transform of the initial condition 9 (x). 1 Assuming 
the spatial mean of 9 vanishes, the variance (6) is 

Var = j 6 2 { Xl t) d 3 x = J \6(k , t)\ 2 d 3 fc , 

which from (13b) becomes 

Var = y*|(?o(fco)| 2 exp |-2k^ (T s • fc ) 2 dsj d 3 fc . (15) 

We thus have a full solution of the advection-diffusion equation for the case of 
a linear velocity field and found the time-evolution of the variance. But what 
can be gleaned from it? We shall look at some special cases in the following 
section. 



2.2 Straining Flow in 2D 

We now take an even more idealised approach: consider the case where the 
velocity gradient matrix a is constant. Furthermore, let us restrict ourselves 
to two-dimensional flows. After a coordinate change, the traceless matrix a 
can only take two possible forms, 

Case (2a) is a purely straining flow that stretches exponentially in one direc- 
tion, and contracts in the other. Case (2b) is a linear shear flow in the x\ 
direction. We assume without loss of generality that A > and U' > 0. The 
form er( 2tl ) is known as the Jordan canonical form, and can only occur for 
degenerate eigenvalues. Since by incompressibility the sum of these identical 
eigenvalues must vanish, they must both vanish. The corresponding funda- 
mental matrices 7 t = exp(— at) are 

- -(_;.«!)■ <»> 

These are easy to compute: in the first instance one merely exponentiates the 
diagonal elements, in the second the exponential power series terminates after 
two terms, because the square of is zero. 



x We are using the convention 



»(fc)=(2^/^)e" a, -"d <l x, 
6(x) = J §{k) e ik x d d k, 



for the Fourier transform in d dimensions. 
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Let us consider Case (2a) , a flow with constant stretching (the case consid- 
ered by Batchelor [1]). The action of the fundamental matrix on fe for Case 
(2a) is 

T t (2a) .fc =(e- At fcoi,e At fco 2 ), (18) 

with norm 

(Tf a) .fc ) 2 = e- 2At fcoi+e 2At fc ^. (19) 

The wavevector k(t) = T^ 2 "' • feo grows exponentially in time, which means 
that the length scale is becoming very small. This only occurs in the direc- 
tion X2, which is sensible because that direction corresponds to a contract- 
ing flow. Picture a curtain being closed: the bunching up of the fabric into 
tight folds is analogous to the contraction. (Of course, it is difficult to close 
a curtain exponentially quickly forever!) The component of the wavevector in 
the x\ direction decreases in magnitude, which corresponds to the opening of 
a curtain. 

Let's see what happens to one Fourier mode. By inserting (19) in (13b), 
we have 

§(ko,t) = ^ (fco)cxp|- K ^ (c- 2As fco 2 +e 2As fc ^)d S | . 

The time-integral can be done explicitly, and we find 

§(k ,t) = o (fe o )cxp{-^ ((e 2At - 1) k 2 2 - (c- 2At - 1) kol)} . 

For moderately long times (t > A -1 ), we can surely neglect e _2A * compared 
to 1, and 1 compared to e 2At , 

0(feo,t)~0o(feo)cxp{-^(e 2At fc o 2 + fc o 2 )} . (20) 

Actually, this assumption of moderately long time is easily justified physically. 
If nk 2 /X <C 1, where k is the largest initial wavenumber (that is, the smallest 
initial scale), then the argument of the exponential in (20) is small, unless 

e 2At > Pe (21) 

where the Peclet number is 

"< = ^- < 22 > 

Thus the assumption that e 2At is large is a consequence of Pe being large, since 
otherwise the exponential in (20) is near unity and can be ignored — variance 
is approximately constant. We can turn (21) into a requirement on the time, 

At^logPc 1 / 2 . (23) 

It is clear from (23) that A -1 sets the time scale for the argument of the expo- 
nential in (20) to become important. The Peclet number influences this time 
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scale only weakly (logarithmically). This is probably the most important phys- 
ical fact about chaotic mixing: Small diffusivity has only a logarithmic effect. 
Thus vigorous stirring always has a chance to overcome a small diffusivity, no 
matter how small: we need just stir a bit longer. 
Note that the variance is given by 

Var = j |0 o (fco)| 2 cxp{-^(c 2At k % + k \) } d 2 k Q , 

which is approximately constant for t <C A _1 logPe 1 ^ 2 . This does not mean 
that the concentration field 

6(x,t) = j 9(k ,t)e [k< - t > x d 2 k (24) 

is constant, even if 6(ko,t) is, because is a function of time 

from (18). This time dependence becomes important for t > A -1 . 

The Peclet number may be thought of as the ratio of the advection time 
of the flow to the diffusion time for the scalar. It is usually written as 

Pe:=UL/n, (25) 

where U is a typical velocity and L a typical length scale. Our velocity estimate 
in (22) is A/fc, and our length scale is k, which are both natural for the problem 
at hand. Just like large Sc, large Pe is the natural setting for chaotic advection. 
In fact, if Pe is small then diffusion is faster than advection, and stirring is not 
really required! Large Pe means that diffusion by itself is not very effective, 
so that stirring is required. We shall always assume that Pe is large. 

We return to (20): the striking thing about that equation is its prediction 
for the rate of decay of the concentration field. Roughly speaking, (20) predicts 

6{x,t) ~ exp{-Pe~V At } (26) 

for At 3> 1. Eq. (26) is the exponential of an exponential — a superexponential 
decay. This is extremely fast decay. In fact, unnaturally so: it is hard to imagine 
a physically sensible system that could mix this quickly. Something more has 
to be at work here. 

If we examine (20) closely, we see that the culprit is the term 

e 2M k t, (27) 

which grows exponentially fast. This term has its origin in the Laplacian 
in the advection-diffusion equation (1): the contracting direction of the flow 
(the xi direction) leads to an exponential increase in the wavenumber via 
the curtain-closing mechanism. This is exactly the mechanism for enhanced 
mixing we advertised on p. 3: very large gradients of concentration are being 
created, exponentially fast. This mechanism is just acting too quickly for our 
taste! 
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So what's the problem? We are doing the wrong thing to obtain our es- 
timate (26). This estimate tells us how fast a typical wavevector decays, and 
it says that this occurs very quickly. What we really want to know is what 
modes survive superexponential decay the longest, and at what rate they de- 
cay. Clearly the concentration in most wavenumbers gets annihilated almost 
instantly, once the condition (23) is satisfied. But a small number remains: 
those are the modes with wavevector closely aligned to the x\ (stretching) 
direction, or equivalently that have a very small projection on the x 2 (con- 
tracting) direction. To overcome the exponential growth in (27), we require 

fc 02 ~ e~ At , (28) 

that is at any given time we need consider only wavenumbers satisfying (28), 
since the concentration in all the others has long since been wiped out by 
diffusion. The consequence is that the fco 2 integral in (24) is dominated by 
these surviving modes. To see this, we blow up the k® 2 integration by making 
the coordinate change k 02 = k 02 c xt in (24), 

/OO />OC ^ ^ 

dfc 01 / dfc 02 9 (k 01 ,k 02 e- xt ) 
-oo J —oo 

x e ife « *exp{-^(fco2 + fcoi)}, (29) 

The decay factor e~ xt has appeared in front. For small diffusivity, we can 
neglect the feoi term in the exponential (it just smooths out the initial con- 
centration field a little). 2 We can then take the inverse Fourier transform 
of 9 a (k ai ,k 02 c- M ), 

flo(fcoi, ^02 c ~ At ) = (2^yi J ®o(x) exp ^—ik 01 xi - (k a2 e~ At x 2 ) dxi dx 2 

and insert this into (29). We then interchange the order of integration: the fcoi 
integral gives a 5- function, and the fco 2 integral gives a Gaussian. The final 
result is 

/oo 
e (c- xt x 1 ,i 2 )G(x 2 -c- xt i 2 ; £)dx 2 , (30) 
-oo 



where 



G{x;i):= ^^c- x2/2f2 (31) 

is a normalised Gaussian distribution with standard deviation £, and we de- 
fined the length 



2 We require the initial condition to be smooth at small scales. Here's why: for 
small k, koi needs to be large to matter in the argument of the exponential. But 
a smooth 9 decays exponentially with fcoi, so there is no variance in these modes 
anyways. 
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I := y/n/X. 

If the initial concentration decays for large \x2 | (as when we have a single blob 
of dye), then (30) can be simplified to 

/>oo 
6» (e~ At x 1 ,i 2 )di 2 . (32) 
-oo 

So the x\ dependence in (32) is given by the "stretched" initial distribution, 
averaged over xi . The important thing to notice is that 

6{x, t) ~ e~ At . (33) 

This is a much more reasonable estimate for the decay of concentration 
than (26)! The concentration thus decays exponentially at a rate given by 
the rate-of-strain (or stretching rate) in our flow. The exponential decay is 
entirely due to the narrowing of the domain for eligible (i.e., nondecayed) 
modes. This "domain of eligibility" is also known as the cone or the cone of 
safety [4, 14]. (In two dimensions it is more properly called a wedge.) The con- 
centration associated with wavevectors that fit within this cone is temporarily 
shielded from being diffusively wiped out, but as the aperture of the cone is 
shrinking exponentially more and more modes leave the safety of the cone as 
time progresses. 

Notice that (33) is independent of k. This brings us to the second most 
important physical fact about chaotic mixing: The asymptotic decay rate of 
the concentration field tends to be independent of diffusivity. But note that a 
nonzero diffusivity is crucial in forcing the alignment (28). The only effect of 
the diffusivity is to lengthen the wait before exponential decay sets in, as given 
by the estimate (23). But this effect is only logarithmic in the diffusivity. 

We can also try to think of (32) in real rather than Fourier space. Con- 
sider an initial distribution of concentration. Our straining flow will stretch 
this distribution in the x\ direction, and contract it in the x 2 direction. Gra- 
dients in x 2 will thus become very large, so that eventually diffusion will limit 
further contraction in the x 2 direction and the distribution will stabilise with 
width \J k/X (see Fig. 1). This is what the Gaussian prefactor in (32) is telling 
us: the asymptotic distribution has "forgotten" its initial shape in x 2 . We say 
that the contracting direction has been stabilised. 

2.3 Shear Dispersion in 2D 

So far we have only considered case (2a) in (16). For case (2b), we have 
from (17) 

T t (2 " ) .fc = (fc 01 ,fc 02 -(7'U 01 ) 

with norm 

(T t (2b) • fc ) 2 = fco? + (fc 02 - U'tk Q1 f . (34) 
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Fig. 1. A patch of dye in a uniform straining flow. The amplitude of the concen- 
tration field decreases exponentially with time. The length of the filament increases 
exponentially, whilst its width is stabilised at I = 're/A. 



Inserting (34) in (14), we have 

9(x,t) = y"^fco)e ifc(t) -*exp (k a \ + (k 02 - U's k Ql f) dsj d 2 k Q . 

We can then explicitly do the time integral in the exponential, 

9(x,t) = 1 0o(k o )e ik W- x 



x exp |-« kol t - ((U't k 0l - k 02 y + k 3 2 ) | d 2 fc . (35) 

The enhancement to diffusion in this case is reflected in the cubic power of 
time in the exponential. This is not as strong as the exponential enhancement 
of case (2a), but is nevertheless very significant. This phenomenon is known 
as shear dispersion or Taylor dispersion. The mechanism is often called the 
Venetian blind effect. Assuming the initial distribution 9q depends only on koi, 
then lines of constant concentration which are initially vertical are tilted by the 
shear flow, in a manner reminiscent of Venetian blinds. The distance between 
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the lines of constant concentration decreases with time as (U't) -1 , which gives 
an effective enhancement to diffusion. The time required to overcome a weak 
diffusivity is thus 

U't > (kolK/U')- 1/3 . (36) 

If we use koi 1 as a length scale and U' as a time scale, we can define a Peclet 
number Pc := V /(fcoi 2 K ) an d rewrite (36) as 

U't > Pe 1/3 (37) 

which should be compared to (23), the corresponding expression for the case 
(2a). Here there is a power law dependence on the Peclet number, rather 
than logarithmic, so we may have to wait a long time for diffusion to become 
important. This makes the linear velocity field approximation more likely to 
break down. 

Let us consider the time-asymptotic limit U't » 1: we might then be 
tempted to neglect everything but the U't fcoi term in the argument of the ex- 
ponential in (35). However, this would be a mistake. To see more clearly what 
happens, define the dimensionless time r := U't and the length x 2 = k/U'. 
Equation (35) then becomes 

9(x,t) = J e (k a )c ik ^- x 

x exp {-x 2 (fc 2 r + fc 2 r + ±fc 2 r 3 - k 01 k 02 T 2 ) } d 2 k . (38) 

The first two terms in the exponential are just what is expected of regular 
diffusion in the absence of flow. The next term is the enhancement to diffusion 
along the x\ direction: it will force the modes fcoi ~ t~ 3 / 2 to be dominant, 
since everything else will be damped away. Similarly, the last term forces fco2 ~ 
r -1 / 2 . Assuming these scalings, the only term that can be neglected for r ^> 1 
is the very first one, fco 2 r. 

We make the change of variable fcoi = T 3 / 2 feoi> fco2 = r 1//2 fc 02 in (38), 

0(X, t) = T- 2 J 9 Q (k Q1 T" 3 / 2 , fc 02 T- 1 ' 2 ) e i(r- 3/2 x 1 -r- 1 / 2 x 2 )fc 01 +ir- 1 / 2 x 2 fc 02 

x exp j-X 2 {ka\ + |^0i - ^01^02) } dfcoi dfc 02 • (39) 

If we approximate #0(^01 t~ 3 / 2 , fc 02 r~ 1//2 ) ~ #o(0, 0), we can do the integrals 
in (39) and find 

9(x,t) * 2V3n X - 2 T- 2 9 (0,0) exp j- ^ ~ 3a ^ T + ^ | . (40) 
For moderate values of x\ (x\ <C x T ), we have 



9(x,t) ~ 2V37rx _2 T- 2 e^ :E ' /x2r o (O,O). 



(41) 
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The width in the x 2 direction of an initial distribution thus increases as x 1 " 1 ^ 2 
\ : ' . This is independent of U' and is exactly the same as expected from pure 
diffusion. The width in the x\ direction in (40) increases as xt^I 2 = U't^fnt 
(see Fig. 2). 



t = o H 6 t=2 




Fig. 2. A patch of dye in a uniform shearing flow. The amplitude of the concen- 
tration field decreases algebraically with time as t~ 2 . The length of the filament 
increases as i 3//2 , whilst its width increases as i 1//2 . 



2.4 Three Dimensions 

In three dimensions, there are three basic forms for the matrix a: 

/Ai \ / 0\ /AO \ 

= A a ; = \ W ; = \ U' A , 

\0 -A1-A2/ \ J7' 0/ \0 0-2A/ 

with corresponding fundamental matrices 
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/e" Al * \ / 0\ 

T< 3a ) = e" A2t ; T< 3b ) = -U't ; (42) 

V c^+^J \\{U'tf -U'tO) 

I e" At \ 
J( 3c ) = -U'te~ Xt e" At . (43) 
V e 2A V 

We can assume without loss of generality that Xi > 0, Ai > A 2 , and 
U' > 0, but the sign of A 2 and A is arbitrary; however, we must have 
A3 = — Ai — A 2 < 0. The case of greatest interest to us is (3a). The rele- 
vant k(t), corresponding to (18), is 

?f a) • fco = (V Alt fcoi , e- A2t fc 02 , e |A3|t fc 03 ) , (44) 
which is used in (14) to give 

9(x,t) = J o (feo)e ife W a; exp{-I K (A- 1 (l-e- 2Alt )fco? 

+ A^ 1 (1 - e" 2A ^) k Q t + IA3I- 1 (c 2 ' A ^ - l) fc 2 ) } d 3 fc . (45) 

What happens next depends on the sign of A 2 : the question is whether e _2A2 * 
grows or decays for t 3> |A 2 | _1 . If A 2 > 0, then we have 



0(x,t)~ J o (feo)e ife( * ) a, exp{-i K (A- 1 fcoi 

+ A 2 - 1 fc 2 + |A 3 r 1 e 2 l A 3l t fc 2 )}d 3 fc , (46) 

whilst for A 2 < 

9(x,t)~ J o (fco)c ife( * ) x cxp{-i K (Ar 1 fcoi 



+ 



M-V^Kfco 2 + IA3I- 1 e 2 l A ^fc 2 ) } d 3 fc . (47) 



Both approximations are valid when t 3> max(A x 1 , |A 2 | x ). For A 2 = the 
situation is similar to the two-dimensional case (2a): 

J 0o(feo)c ifc(t) K cxp{-^-(fc o 2 +c 2Alt fc o 2 )}d 3 fc o , 

valid when t ^> A^f 1 . 

The rest of the calculation is very similar to the two-dimensional case (2a), 
in going from (24) to (32). In both (46) and (47) the X3 direction is stabilised, 
that is we need to blow up the fco3 integral to remove the time dependence 
from the exponential, and find that the integral is dominated by ko 3 — 0. 
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The x 2 direction is also stabilised in (47), so we can set k 02 — 0. We thus find 
for A 2 > 0, 

0(x, t) ~ e-' A3 l* G{x z ;l 3 ) J 9 a ( e - Xlt x 1 ,e- X2t x 2 ,x 3 ) dx 3 , (48) 
and for A2 < 0, 

9(x,t) ~ c-^+^t G(x 2 ;e 2 ) G(x 3 ;l 3 ) J O^^xx, x 2 , £3) dx 2 dx 3 , 

(49) 

where li :— ^ k/W\. Contracting directions have their spatial dependence 
given by a time-independent Gaussian, with an overall exponential decay; 
stretching directions do just that: they stretch the initial distribution, with 
no diffusive effect. Solutions of the form (48) are called pancakes, and those 
of the form (49) are called ropes or tubes. 

There is another way of thinking about the asymptotic forms (32), (48), 
and (49) [15]: contracting directions are stabilised near some constant width 
tj, and expanding directions lead to exponential growth of the width of an 
initial distribution along the direction. Thus, the volume of the initial distri- 
bution grows exponentially at a rate given by the sum of A,'s associated with 
stretching directions, but the total amount of 9 remains fixed (the mean is 
conserved). Hence, the concentration at a point should decay inversely pro- 
portional to the volume, which is exactly what (32), (48), and (49) predict. 



3 Random Strain Models 

In Section 2 we analysed the deformation of a patch of concentration field (a 
'blob') in a linear velocity field. Though this is interesting in itself, it is a far 
cry from reality. We will now inch slightly closer to the real world by giving a 
random time dependence to our velocity field. 

3.1 A Single Blob 

Consider a single blob in a two-dimensional linear velocity field of the type we 
treated in Section 2.2 (case (2a)). Now assume the orientation and stretching 
rate A of the straining flow change randomly every time r. This situation is 
depicted schematically in Figure 3. We assume that the time r is much larger 
than a typical stretching rate A, so that there is sufficient time for the blob to 
be deformed into its asymptotic form (32) at each period, which predicts that 
at each period the concentration field will decrease by a factor exp(— A^r), 
where A^ is the stretching rate at the ith period. The concentration field 
after n periods will thus be proportional to the product of decay factors, 

0~e-* ( %- AC2>T --.e- A(n \ 

= c -(A<i> + A<2> + ...+A<"> K (5Q) 
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e -AWr e -»«r e -A<3). 




Fig. 3. A single blob being stretched for a time r by successive random straining 
flows. The amplitude of the concentration field decays by exp(— A^V) at each period. 

We may rewrite this as 

9~e~ Ant , (51) 

where t — nr, and 

1 11 
n t— i 

is the 'running' mean value of the stretching rate at the nth period. As we let n 
become large, how de we expect the concentration field to decay? We might 
expect that it would decay at the mean value A of the stretching rates A^ . 
This is not the case: the running mean (51) does not converge to the mean A. 
Rather, by the central limit theorem its expected value is A, but its fluctuations 
around that value are proportional to 1/ y/i . These fluctuations have an impact 
on the decay rate of 8. 

The ensemble of variables A^ is known as a realisation. Now let us imagine 
performing our blob experiment several times, and averaging the resulting 
concentration fields: this is known as an ensemble average over realisations. 
Ensemble-averaging smooths out fluctuations present in each given realisation. 
We may then replace the running mean A n by a sample-space variable A, 
together with its probability distribution P(A, t). The mean (expected value) 
of the ath power of the concentration field is then proportional to 

/•oo 

/ e- aAt P(A,t)dA. (52) 
Jo 
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The overbar denotes the expected value. The factor c~ aAt gives the amplitude 
of 6 a given that the mean stretching rate at time t is A, and P(A, t) measures 
the probability of that value of A occuring at time t. 

The form of the probability distribution function (PDF) P(A, t) is given 
by the central limit theorem: 

P(A,t)~G(A-A;^pTt), (53) 

that is, a Gaussian distribution (31) with mean A and standard devia- 
tion WP/t. The quantity /3 is a measure of the nonuniformity (or fluctua- 
tions) of stretching in the flow, both spatially and temporally. The decrease 
of the standard deviation \J (3/t with time reflects the convergence of the av- 
erage stretching rate to the Lyapunov exponent A. I will say more about (3 in 
Section 4.1, when we look at a practical example. 

Actually, the central limit theorem only applies to values of A that do not 
deviate too much from the mean. The theorem understimates the probability 
of rare events; a more general form of the PDF of A comes from large deviation 
theory [16,17], 

P(i )t )~^M e - t8 ("). (54) 

(A derivation of (54) is given in Appendix C.) The function S(x) is known as 
the rate function, the entropy function, or the Cramer function, depending on 
the context (that is, which literature one is reading). It is a time- independent 
convex function with a minimum value of at 0: 8(0) = §'(0) = 0. If A is near 
the mean, we have 

B(A-A) ~ ±S"(0)(/L-i) 2 , (55) 

which recovers the Gaussian result (53) with (3 — 1/S"(0). Both (53) and (54) 
are only valid for large t (which in our case means t » r, or equivalently n ^> 1). 
We can now evaluate the integral (52) with the PDF (54), 

/>OC 

W~ / c - tH ^ dA-e-^, (56) 
Jo 

where we have omitted the nonexponential prefactors, and defined 

H(A) := aA + §(A-A). 

Since t is large, the integral is dominated by the minimum value of H{A): 
this is the perfect setting for the well-known saddle-point approximation. The 
minimum occurs at A sp where H'(A sp ) = a + §'(A sp — A) = 0, and is unique 
because § is convex and has a unique minimum. The decay rate is then given 
by 

la =H(A sp ), with H'(A sp )=0. (57) 

There's a caveat to this: for a large enough the saddle point A sp is nega- 
tive. This is not possible: the stretching rates are defined to be nonnegative 
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(the integral (56) involves only nonnegative A). Hence, the best we can do 
is to choose A sp = — the integral (56) is dominated by realisations with no 
stretching. Thus, in that case -f a = H(0), or 

7 « = §{-A). (58) 

We re-emphasise: for small enough a, the saddle point is positive and the 
decay rate is given by (57). Beyond that, we must choose zero as the saddle 
point and the decay rate is given by (58) . To find the critical value a cr ; t where 
we pass from (57) to (58), observe that this happens as the saddle point nears 
zero. Thus, we may solve our saddle point equation H'(A sp ) = by Taylor 
expansion, 

H'(A sp ) ~ a crit + S'(-A) + A sp §"(-A) = 0. (59) 

But the saddle point will not be small unless the first terms cancel in (59), 
that is Qfcrit = —§'(—A). We may thus recapitulate the result for the decay 
rate, 

(aA sp + §(A sp -A), a<-S'{-A); 

7a = { _ ' (60) 

[§(-A), a>-S'(-A). 

Clearly j a is continous, and it can be easily shown that d-y a /da is also con- 
tinuous. 

As an illustration, we use the Gaussian approximation (55) for the Cramer 
function, with = 1/S"(0). The critical a is a crit = -S'(--A) = A/0. The 
saddle point is positive for a < A/0, so from (60) we get 

a(A-\a0), a<A/0; 

-I - ( 61 ) 

A 2 /20, a > A/0. 

This is plotted in Figure 4. Notice that the solid curve (for a random flow) 



3 




-2 -1 1 2 3 4 5 6 



a 

Fig. 4. Decay rate (61) for the concentration of a blob in a Gaussian random 
stretching flow (solid curve). The dashed line is for a fixed, nonrandom flow as in 
Section 2.2. Here A = 1, = 1/4, so a cr it = A/0 = 4. 



20 Jean-Luc Thiffeault 



lies below the dashed line (for a nonrandom flow). This is a general result: 
if f(x) is a convex function and x a random variable, Jensen's inequality says 
that 

W)>f{x). (62) 
Now, e~ atA is a convex funtion of A, so we have 

i 3 ^ 5 " > e"" a , 

which means that the rate of decay satisfies 

7a < a A 

which is exactly what is seen in Figure 4. Thus, fluctuations in A inevitably 
lead to a slower decay rate "f a . 

Stronger fluctuations also means that the decay rate 7 Q saturates more 
quickly with a. Clearly, in the absence of fluctuations we recover the nonran- 
dom result: A/p is infinite and only the a < A/ f3 case is needed in Eq. (61). 
If there are lots of fluctuations, A//3 is small, and there is a greater probabil- 
ity of obtianing a realisation with no stretching. For large enough fluctuations 
this exponentially-decreasing probability dominates, and we obtain the second 
case in (60). 

3.2 Many Blobs 

In Section 3.1 we considered the evolution of the concentration of a single 
blob of concentration in a random straining field. Now we turn our atten- 
tion to a large number of blobs, homogeneously and isotropically distributed, 
with random concentrations. We assume that the mean concentration over 
all the blobs is zero. A simplified view of this initial situation is depicted in 
Figure 5(a), with shades of gray indicating different concentrations. If we now 
apply a uniform straining flow of the type (2a) (see Section 2.2), the blobs 
are all stretched horizontally (the x\ direction) and contracted in the vertical 
{x-i) direction, as shown in Figure 5(b). They are pressed together in the x^ 
direction until diffusion becomes important (Figure 5(c)). The effect of diffu- 
sion is to homogenise the concentration field until it reaches a value which is 
the average of the concentration of the individual blobs. This is depicted by 
the long gray blob in Figure 5(d), which will itself keep contracting until it 
reaches the diffusive length I. 

Of course, here the initial concentration field 9q represents the concentra- 
tion of all the homogeneously-distributed blobs together, so it does not decay 
at infinity: we must thus use Eq. (30) rather than (32). The summing (and 
hence averaging) over blobs is manifest in Eq. (30), which contains an integral 
over the initial distribution 8q in the x-i direction, windowed by a Gaussian. 

In practice, this implies that the expected value of the concentration at a 
point x on the gray filament is given by 
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(a) (b) 




Fig. 5. (a) An initial distribution of blobs with random concentrations (b) are 
stretched by a constant strain (c) until they reach the diffusive limit in the contract- 
ing direction and begin to overlap, (d) Finally, they combine into one very long blob 
with the average concentration of all the blobs. 

(0(M))biobs~e-^^^O, (63) 

i 

where 9$ is the initial concentration of the ith blob, and (-) blobs denotes 
the expected value of the sum over the overlapping blobs at point x (not the 
same as spatial integration (•)). We assume that N ^> 1 blobs have overlapped. 
Equation (63) gives the concentration at a point, summed over TV overlapping 
blobs. Of course, Eq. (63) converges to zero for large N, because the blobs 
average out. Not so for the fluctuations at that point: by the central limit 
theorem, we have 



<*»(*, <)) blobs ~ a" 9 * £ 6f = 7Ve _2ilt 61 , (64) 

i 

since the initial blobs have identical distributions. The blob-summed fluctua- 
tion amplitude (# 2 ) blobs is thus proportional to the number N of overlapping 
blobs. But the number of overlapping blobs is proportional to e At : as time 
increases more and more blobs converge to a given x in the contracting di- 
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rection and overlap diffusively (this can be seen in Eq. (30): the width of the 
windowing region grows as e At ). Assuming the variance of 9$ is finite, we 
conclude from (64) that 

(^))HL~*- At/2 - (65) 

Compare this to (51) for the single-blob case: the overlap between blobs has 
led to an extra square root. Thus, the ensemble averages (6 2 (x, t))^ lohs for the 
overlapping blobs are computed exactly as in Section 3.1, resulting in (60). 
Because of the assumption of homogeneity, the point-average is the same as 
the average over the whole domain (see Section 4 for more on this), and we 

have 3 

<0TH0 2 M> blobs ~e-^ 4 , (66) 

with j a given by (60). (In (66) the angle brackets denote spatial averaging, 
not spatial integration, because the total variance is infinite in this case.) 

3.3 Three Dimensions 

In three dimensions, we will only treat Case (3a) (a purely straining flow) of 
Section 2.4. For A 2 < 0, where the asymptotic concentration is given by (49) 
(ropes), the situation is basically identical to the 2D case of Sections 3.1- 
3.2: the statistics of the stretching direction Ai determine j a from (60). The 
contracting directions x-i and xz are stabilised by diffusion. 

For A 2 > 0, the asymptotic concentration is given by (48) (pancakes). 
We have two fluctuating quantities to worry about (Ai and A 2 ). But since the 
decay rate in (48) only depends on A3 , we can instead focus on its fluctuations 
only. For a single blob, the average (52) is then replaced by 

/>oo 

/ e-^l^d^l.tJdlAa^e-T-*, (67) 
Jo 

where of course A 3 is the average of A3. This PDF achieves a distribution of the 
large-deviation form (54). The analysis thus follows exactly as in Sections 3.1 
and 3.2, except the Cramer function for |yl 3 | must be used. 4 

4 Practical Considerations 

One may rightly wonder if the blobs in a random uniform straining flow de- 
picted in Section 3 bear any resemblance to reality. The single-blob scenario 

3 In going from (65) to (66), we've implicitly assumed that the initial concentra- 
tion field has Gaussian statistics, because we've used the fact that the higher even 
moments are proportional to powers of the second moment. 

4 There are a few exceptional cases to consider [15]. 
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doesn't, but the many-blobs scenario has a fighting chance, as we will try to 
justify here. There are two important considerations: where docs the ensemble- 
averaging come from, and what are the stretching rates given by? 

The decay rate (60) depends crucially on ensemble-averaging: with that 
averaging the decay rate fluctuates wildly for a given realisation. At the end 
of Section 3.2 we assumed that homogeneity allowed us to generalise from the 
average at a point to the average over the whole domain. But the average over 
the whole domain can actually do a lot more for us: it can provide the ensemble 
of blobs that we need for averaging! Thus, we can forget about speaking of 
realisations as if we were running many parallel experiments, and instead 
speak of the moments of the concentration field as given by an average over 
randomly-distributed blobs. The decay rate will then be naturally smoothed- 
out over blobs experiencing different stretching histories. The saturation of the 
decay rate with a in Eq. (60) is due to 9 2a being dominated by the fraction 
of blobs that have experienced no stretching. 

What about the stretching rates A? Luckily, it is not them but their time- 
average A that matters. If we imagine following a blob as it moves through 
the flow, we can see that this time-averaged stretching rate is nothing but the 
finite-time Lyapunov exponent associated with this blob and its particular 
initial condition. A given blob will be constantly reoriented as it moves along 
in the flow, so its finite-time Lyapunov exponent is not just the average of the 
stretching rates (in fact, it must be strictly less than this average). But in a 
chaotic system we are guaranteed that, on average, these reorientations do not 
lead to a vanishing (infinite-time) Lyapunov exponent. This is guaranteed by 
the celebrated Oseledec multiplicative theorem for random matrices [18]. 5 We 
may thus use for P(A,t) the distribution of finite-time Lyapunov exponents, 
which is well-known to have the large-deviation form (54) [19]. 

The result of these considerations is the local theory of passive scalar de- 
cay It is called local because of the reliance of such a local concept as the 
finite-time Lyapunov exponents, which come from a linearisation near fluid 
element trajectories. In Section 4.1 we discuss a specific example. We post- 
pone a discussion of the validity of the local theory until Section 4.2, but for 
now we point out that it is known to be exact at least in some simple model 
flows [15,20]. 

The derivation presented in this Section was based on the work of Balkovsky 
and Fouxon [15], who used a slightly more rigorous approach. Son [21] also 
obtained the decay rate (60) using path-integral methods. Earlier, Anton- 
sen et al. [8] derived the decay rate for the second moment (# 2 ) in terms of 
the Cramer function, using a different (and not quite equivalent) approach, 
though they did not allow for the second case in (60). 



5 The reorientations also tend to decrease the correlation time r [15]. 
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4.1 An Example: Flow in a MicroChannel 

We illustrate how to compute the decay rates j a with a practical problem. 
Specifically, we will use a three-dimensional model of a microchannel. The sys- 
tem is shown in Figure 6. It consists of a narrow channel, roughly 100 wide 




Fig. 6. MicroChannel with a periodic patterned electro-osmotic potential at the 
bottom. The arrows indicate the direction of fluid motion at the bottom. The width 
of the channel is about 100 ^m and its height 10-50 /J,m, and the period of the pattern 
is L. A typical mean fluid velocity is 10 2 -10 3 /wn/s. 

and slightly shallower. These types of channel are widely used in microfluidics 
applications ( "lab-on-a-chip" ) , and often one wants to achieve good mixing 
in the lateral cross-section of the channel. This is difficult, since the Reynolds 
number of the flow varies between 0.1 and 100 — far from turbulent. Clever 
techniques have to be used to induce chaotic motion of the fluid particle tra- 
jectories in order to enhance mixing. Stroock et al. [22] used patterned grooves 
at the bottom of the channel to induce vortical motions, and found that the 
mixing efficiency was dramatically increased. Here we use a variation on this 
where the bottom is pattern with an electro-osmotic coating, which induces 
fluid motion near the wall [23] . The effect of the electro-osmotic coating is well- 
approximated by a moving wall boundary condition. The pattern is chosen 
in a so-called herringbone pattern to maximise the mixing efficiency (though 
not in a staggered herringbone, which is even better but is more difficult to 
model). Rather than solving the full equations numerically, we adopt here an 
analytical model based on Stokes flow in a shallow layer [24]. The longitudi- 
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rial (x) direction is taken to be periodic. The flow is steady, but because it is 
three-dimensional it can still exhibit chaos. 

Figure 7 shows two Poincare sections for the flow. These are taken at two 
constant x planes, one at x = and the other at the mipoint of the a;-periodic 
pattern. The two colours represent two trajectories that have periodically 
punctured those planes many times over. It is clear from the Figures that the 
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Fig. 7. Poincare sections for the microchannel. The red and blue dots represent the 
same trajectory periodically puncturing two vertical planes many times over (blue 
if in the same direction as the flow, red otherwise). The green and yellow dots show 
two trajectories in regular, nonmixing regions. 



flow contains large chaotic regions, as well as smaller regular regions (known 
as islands). We focus here on the chaotic regions. 

Now that we have established (or at least strongly suspect) the existence 
of chaotic regions, we can compute the distribution of finite-time Lyapunov 
exponents. There are many ways of doing this: because we are not interested 
in exetremely long times, the most direct route may be used. We have an 
analytical form for the velocity field, so the velocity gradient matrix is eas- 
ily computed. This allows us to linearise about trajectories in the standard 
manner [19,25]. Each trajectory will thus have a finite-time Lyapunov ex- 
ponent associated with it, which shows the tendency of infinitesimally close 
trajectories to diverge exponentially. This is then repeated over many different 
trajectories within the same chaotic region, and a histogram is made of the 
finite-time Lyapunov exponents. This histogram changes with time, as shown 
in Figure 8. For these relatively early times, it changes dramatically and does 
not show a self-similar form. 

The evolution of the mean and standard deviation of the distribution is 
shown in Figure 9. The mean is converging to a constant A ~ 0.116, and the 
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Fig. 8. Evolution of the distribution of finite-time Lyapunov exponents for the 
microchannel. The average crossing time for particles in the channel is L/U. 

standard deviation is decreasing as W '/3/t, with (3 ~ 0.168. (Both A and (3 are 




Fig. 9. (a) Evolution of the mean A of the distribution of Lyapunov exponents. 
The mean converges to A ~ 0.116 s" 1 . (b) Standard deviation of the distribution 
of Lyapunov exponents versurs l/\/t. The straight line represents yj /3/t, with fi ~ 
0.168 s" 1 . 
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fitted values, and arise from the complicated nonlinear nature of fluid parti- 
cle trajectories — they cannot be predicted or calculated from first principles 
except in the simplest cases.) These facts taken together are strongly indica- 
tive that the distribution is converging to a Gaussian of the form (53). This 
is easily confirmed by plotting the PDFs at different times and rescaling the 
horizontal axis by y/i, as shown in Figure 10. Note that this case exhibits 



t=200 




-1 -0.5 0.5 1 



Vi(A-A) 

Fig. 10. Rescaled distribution of finite-time Lyapunov exponents at different times. 
The dashed line is the Gaussian form (53), with parameters as in the caption for 
Figure 9. 

a particularly nice Gaussian form, which is not necessarily the norm for all 
chaotic flows. 

Using the values for A and j3 we just obtained, we can calculate the decay 
rates j a with the Gaussian approximation. The ratio A/ (3 is 0.69, so the change 
in character in (61) occurs at a > 0.69. Since from (66) the decay of (0 2 ^ is 
given by j a , this means that moments of order 2a > 1.38 will decay at the 
same rate. This includes the variance (0 2 ), so we have from (66) 

(6» 2 )~e" 71 *, with 7! = A 2 /2P~ 0.040 s" 1 . 

The mixing time is thus 7-f 1 ~ 25 seconds. This is about a factor of four 
improvement over the purely diffusive time for, say, DNA molecules (k ~ 
10~ 10 m 2 s ). This is not spectacular, but can be greatly increased by stag- 
gering the herringbone pattern. The mixing time assuming the decay proceeds 
at the rate of the mean Lyapunov exponent A is roughly 9 seconds, so that 
the fluctuations multiply this by a factor of three! 
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Of course, we do not know if this is actually a good estimate for the 
mixing time, since we haven't directly solved the advection-diffusion equation 
numerically: this is prohibitive in a three-dimensional domain for such a small 
diffusivity. This is one of the advantages of the local theory: it is usually less 
expensive to compute the distribution of finite-time Lyapunov exponents than 
it is to solve the advection-diffusion equation equation directly. We will say 
more on the validity of the local theory in Section 4.2. 

4.2 Limitations of the Local Theory 

So is this local theory of mixing correct? Well, certainly not always, even in the 
Batchelor regime. There are many assumptions underlying the model, some of 
them difficult to verify (Do blobs really undergo a series of stretching events 
as described here? Do correlations between these events matter?) My feeling 
is that sometimes it will, but most of the time it won't. More experiments 
and numerical simulations are needed to get to the bottom of this. For a 
detailed discussion of possible problems with the local theory, see Fereday 
and Haynes [26] . They make a good case that the theory must break down for 
long times: the blobs discussed here meet the boundaries of the fluid domain 
and must begin to fold. The folding forces them to interact with themselves 
in a correlated fashion. We enter the regime of the strange eigenmode [27], 
which has received a lot of attention lately [14,26,28-35]. Maybe we'll hear 
more about that in ten years. . . . 6 
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A The Advection— Diffusion Equation in a Comoving 
Frame 

We start from the advection-diffusion equation (1) and derive its form (8) for 
a linearised velocity field. We want to transform from the fixed spatial coor- 
dinates x to coordinates r measured from a reference fluid trajectory x (t). 



6 Recently, Tsang et al. [36] have tested the local theories to an astonishing preci- 
sion for the Zeldovich sine flow, whilst Vanneste and Haynes [37] have convincingly 
demonstrated that the local theory holds when the system is dominated by the rem- 
nants of the continuous spectrum of the advection operator, whereas global aspects 
must be considered when the slowest-decaying mode is regular. 
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The coordinates r are not quite material (Lagrangian) coordinates, since we 
follow the trajectory of only one fluid element. 
We thus let 



x = xo(t) + r, 
and write the concentration field as 



dx Q (t) 
dt 



v(x (t),t), 



(68) 



9(x,t) = 6(r,t). 
The time derivative of 9 can be written 



d_ 

dt 



9(x,t) 



d_ 

dt 



9{r,t)+V r 9 



~ dr 

at 



(69) 



where d/dt\ x denotes a derivative with x held constant. Now from (68) 



dr 
8t 



dx . . , . 

-— = -v(xo(t),t). 



Spatial derivatives are unchanged by (68): \7 X 9 = V r 9. Hence, inserting (69) 
into (1), we find 



a_ 

dt 



9 + {v(x (t) +r,t)- v(x (t),t)} -V r 9 = K V 2 J . 



We Taylor expand the velocity field in (70) to get 



d_ 

dt 



r ■ <r(t) ■ V r 9 = nV 2 r 9, <r(t) := Vv(x (t),t), 



(70) 



(71) 



where we neglected terms of order \r\ 2 . This is only valid if the velocity field 
changes little over the region we consider (i.e., if it is smooth enough), which 
is true for large Schmidt number. Equation (71) is the same as (8), and tells 
us how to find a(t). 



B Volume Preservation 

It is useful to know how to prove that a divergence-free vector field will lead 
to volume-preservation of a small blob of integrated trajectories. Mathemat- 
ically, we want to go from (11) to (12). We start from the definition of the 
determinant of a <i-dimensional matrix, 

dctT= — ti x ...i d eji-.-jd Tiiji • • • ^idid (72) 

if -id if jd 
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where we have dropped the T subscript from T for this section, and e is the 
fully-antisymmetric Levi-Civita symbol. Taking a time derivative of (72), 

<9 t dctT = ^ii-i <1 ^i-j li dt7i ljl ---7i djd (73) 

*i ---id h—3i 

since the d terms obtained after using the product rule for derivatives are the 
same. We use the equation of motion (11) for T, 

<9 t dctT=— /j j e ii---id e h---3d ^ht fiji ^2J2 ' ' ' fidjd • 

3l---3d I 

We rearrange the sums on the right, 

<9 t detT= — ^ &i 1 e7gj 1 I ^ ^ ^ e ii--4 d e h—3d ^iih ' ' ' ^idid 

ii,3i,t \i2---id 32— 3d 

and recognise that, up to a factor of detT, the term in the parentheses is the 
cofactor representation of the inverse of T: 

<9 t detT=- ^ <n 1 e7e jl {7" 1 ) iin detT = -^a^eS^e detT 

iijii ii,£ 

so that finally 

<9 t dctT = -(Trcr) detT. 

We conclude that if Trcr = then the determinant of T is constant. Since its 
initial condition is unity (from Eq. (11)), then it must remain so for all time. 



C Large Deviation Theory 

In this Appendix we will justify the large-deviation form of the PDF, Eq. (54), 
assuming little prior knowledge of probability theory. 

First, we define the generating function e~ s ( fc ) of a random variable x by 

e- s(k) = J p(x)c- ikx dx, 

that is, the generating function is simply the Fourier transform of the PDF 
of x. We have s(0) = 0, x = — is'(0), and x 2 — x 2 = s"(0). Now define the 
random variable X to be the mean of n variables, 
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where the independent and identically distributed with PDF p(xi) — 

p(x). How do we find the PDF P{X n ) of X n , in the limit where n is large? 
First observe that (from here on we drop the subscript on X n ) 

P{X) = J p( Xl ) ■ ■ -p(x n ) S ( ^ 1 + " n +Xn - dan ■ ■ ■ dx n 

since the joint PDF p(x\, . . . , x n ) — p{x\) ■ ■ ■ p(x n ) by independence of the a;,. 
The generating function e^ s ^ for P(X) is then 

e s(k) = J P (x)c-' lkX dX 

= j p{x 1 )---p{x n )5( ^ l + " n +Xn -X^jc- ikX dzi • • • dx n dX 

We do the X integral, and then observe that we get a product of n identical 
Xi integrals, each of which is equal to e~ s ( fe /™': 



e S(k) 



= J P (xi) ■ ■ -p{x n ) e -iWn)(x 1 + -+ Xn ) ^ . . . = e -na(fc/n) _ 



Thus the generating function for X is the nth power of the generating function 
for x. We can invert the Fourier transform to find the PDF P(X): 

P(X) = — [ c- s « e ikx dk = — / e -»MK)-iKX) dR (?4) 
2tt J 2n J 

where K — k/n. We let 

H(K,X) = s(K) -iKX; 

For large values of n the integral in (74) is dominated by the stationary 
points K sp (X) of H(K,X) (saddle-point approximation): 



ri M 

K sp (X) such that — {K sp ,X) = s'(K sp )-iX = 0. (75) 



(I will often leave out the X dependence of K sp (X) to shorten the expressions.) 
In that case we can approximate the integrand in (74) using 

H{K,X) = H(K sp , X) + i s"(K sp )(K - K sp f +0({K- K sp f) 

which allows us to do the integral explicitly: 



p(X) = / - -nH(K sp (X),X) (76) 

where K sp (X) is given by (75). As a final step, let us calculate the mean of 
X using this PDF: 
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X = /xP(X)dX = /x ^=g= e-^'^dX. (77) 

Again, for large n we can use the saddle-point method to evaluate this inte- 
gral. The important observation is that the saddle-point X of H sp (K(X), X) 
satisfies 

^(K sp (X ),X ) = ^L(K SP (X ),X ) ^(X ) -iK sp (X ) = 0. 

The dH/dK term vanishes because it is evaluated at K sp ; hence, K sp (Xo) = 0, 
which implies H (K sp (Xa), Xo) = 0. Inserting this into the integral (77), we 
find X = X n : the mean of X and the minimum of H coincide. This means 
that it makes sense to define 

$(X - X) := H(K sp (X), X), with S(0) = and S'(0) = 0, (78) 

which is the sought-after Cramer function. Note also that §"(X — X) = 
1/ s"(K sp (X)) 1 and that for large n the non-exponential coefficient in (76) 
can thus be approximated by evaluating it at the saddle-point K sp (X) = 0, 
with s"(0) = 1/S"(0). The final form of our large-deviation result is thus 



P(I) = ,i!!l c -»s(«) (79) 

V 2ir 

which is the same as Eq. (54). 

As a simple example (treated in every textbook, see for example [17]), 
consider a random variable x with PDF 

p(x) = (1 -s)S(x-x+) +e6(x-x-), (80) 

where x + > x_ are constants — this is a binomial distribution (or Bernoulli 
distribution in this case). If we take the mean X of n such variables, what is 
the PDF of X for large n? First, we compute the generating function for x, 

c- s W = J {(l-e)6(x~x + )+e5(x~x-)}c- [kx dx 

= (l-e)c- ikx + +ec-' lkx - . (81) 

We take the logarithm to obtain s(k) and find K sp (X) by solving the saddle- 
point equation (75), 

dH , „ 1 (1-ex + -X 



dK 



3 '(K sp )-\X = K sp (X) = ^\og(^ 



where A := x + — x_ and we restrict a;_ < X < x + . Inserting this into 
H(K sp (X),X), we find from (78) 

o/v X — x- , (1 — £X + — X\ ( x + - X s 

§(X -X) = log -± + log ' 



X-x-J ° eA 
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It is easy to verify that, since X = (1 — e)x + + ex_, we have 8(0) = §'(0) = 
andX^-X 2 = 1/S"(0) = e(l-e)A 2 . 

The binomial distribution (80) is a useful model of stretching of an in- 
finitesimal line segment by a uniform incompressible straining flow in two di- 
mensions, assuming the straining axes of the flow change direction randomly 
at regular intervals r. If we set x± = ±At = ±/3, where A is the strain rate, 
then X is the averaged logarithm of the length £ of the segment, i.e. £ = e nX . 
Thus, the mth power of the length of the segment will on average grow as 

|^ = ^X =e -S(imn) =e -ns(im) = j ^ _ £ ) + £ e -/3mj" _ (g 2 ) 

We know that £~ 2 must be constant in a 2D incompressible flow [13], so that 
the term in braces in (82) must be unity. We use this to solve for e, 

which then allows us to use (82) to write the growth rate \ m of line segments 
as 

l,-^— 1, / cosIi(to + 1)/3 n 



Xra = l0g^ m =-l0g 

in r \ cosh p 

The Lyapunov exponent, which is given by dxm/dm at m = 0, has a value 
of A tanh /3 for this flow: it is less than for a uniform straining flow because of 
the time taken for the segment to realign with the new straining axis when 
its direction changes. 
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